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Abstract. Evolutionary solar models are computed using a 
new stellar evolution code, MoSEC (Modular Stellar Evolution 
Code). This code has been designed with carefully controlled 
truncation errors in order to achieve a precision which reflects 
the increasingly accurate determination of solar interior struc- 
ture by helioseismology. A series of models is constructed to 
investigate the effects of the choice of equation of state (OPAL 
or MHD-E, the latter being a version of the MHD equation 
of state recalculated by the author), the inclusion of helium and 
heavy-element settling and diffusion, and the inclusion of a sim- 
ple model of mixing associated with the solar tachocline. The 
neutrino flux predictions are discussed, while the sound speed 
of the computed models is compared to that of the sun via the 
latest inversion of SOI-MDI p-mode frequency data. The com- 
parison between models calculated with the OPAL and MHD-E 
equations of state is particularly interesting because the MHD-E 
equation of state includes relativistic effects for the electrons, 
whereas neither MHD nor OPAL do. This has a significant ef- 
fect on the sound speed of the computed model, worsening the 
agreement with the solar sound speed. Using the OPAL equa- 
tion of state and including the settling and diffusion of helium 
and heavy elements produces agreement in sound speed with 
the helioseismic results to within about ±0.2%; the inclusion 
of mixing slightly improves the agreement. 
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1. Introduction 

Most of the differences between observed p-mode frequencies 
and those calculated from evolutionary solar models are be- 
lieved to reflect uncertainties in the ingredients of solar mod- 
elling. Much effort has consequently been applied to finding 
which aspects of the uncertain physical assumptions are respon- 
sible for the differences between observation and theory (e.g. 
Turck-Chieze & Lopes 1993). This has led, for example, to the 
use of the MHD equation of state (Mihalas et al. 1988), the 
improved calculation of opacity in the OPAL tables (Iglesias 
& Rogers 1991; Rogers & Iglesias 1992), the inclusion of the 


effects of helium settling and diffusion (Bahcall & Pinsonneault 
1992b; Christensen-Dalsgaard et al. 1993), and finally the de- 
velopment of models with heavy element settling and diffusion 
(Proffitt 1994; Bahcall & Pinsonneault 1995). 

Several components go into a typical recipe for solar evo- 
lution: these include a set of opacity tables, a set of equation 
of state tables, and a set of reaction rates and temperature de- 
pendences for each appropriate nuclear reaction. Additional as- 
sumptions concern the choice of convection theory used, the 
model of the solar atmosphere employed, and various input pa- 
rameters, such as present-day solar age and luminosity. Given 
a certain set of choices for each of these ingredients, any evolu- 
tion calculation should ideally give the same result to within the 
numerical accuracy of its integration scheme. Much progress 
has been made towards attaining this goal (e.g. Christensen- 
Dalsgaard 1991a) by eliminating the errors in complex evolu- 
tion codes which previously caused them to disagree. 

Given the accuracy with which p-mode frequencies are cur- 
rently being measured by such instruments as SOI-MDI (Koso- 
vichev et al. 1997) and GONG (Harvey et al. 1996), subtle as- 
pects of the input physics are becoming subject to scrutiny. In 
order to facilitate such progress, not only should errors be elim- 
inated in evolution codes, but also codes should be developed 
with very carefully controlled truncation errors, such as the CE- 
SAM code described by Berthomieu et al. (1993). With these 
considerations in mind, we have developed a new stellar evolu- 
tion package named MoSEC, designed with minimization and 
control of truncation errors very' much to the fore. 

MoSEC is used to generate a series of seven evolutionary 
solar models with varying input physics. These are designed 
to investigate the effect of including settling and diffusion of 
helium and heavy elements, the effect of the choice of equa- 
tion of state, the effect of ignoring the variation of the heavy- 
element abundance in the calculation of the equation of state 
(as has been done by various authors), and the effect of mix- 
ing below the base of the convection zone. The neutrino flux 
predictions are compared to measured solar values, while the 
sound-speed is compared to that of the sun using an inversion 
of the latest SOI-MDI p-mode frequencies. In the case of the 
equation of state, MHD-E is a new MHD-like calculation taking 
into account the relativistic correction to the electron pressure, 
a correction which was not included in either the OPAL or the 
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original MHD equations of state. A solar model is calculated 
usino MHD-E equation of state, and compared to two models, 
one calculated with OPAL, and the other also calculated with 
OPAL, but adding in the relativistic correction to the electron 
pressure The last model enables us to evaluate what fraction 
of the difference between the OPAL and MHD-E models arises 
from the different treatment of the electron pressure. 

2. Formulation of the problem 

2.1. General principles 

The term “standard solar model" defines a certain set of simpli- 
fying assumptions under which the calculation of stellar evolu- 
tion may more easily be performed. A uniform initial chemi- 
cal composition is assumed, rotation, magnetic fields and mass 
loss are neglected, and spherically-symmetric hydrostatic equi- 
librium is imposed. Mass is chosen as the independent variable 
rather than radius, since the equations of stellar structure be- 

come less nonlinear. . 

The partial differential equations of stellar evolution are as 
usual expressed as two sets of ordinary differential equations: 
one that describes the equilibrium structure at a particular time 
given the chemical composition, and another that describes the 
evolution with time of the chemical composition. The solution 
of the former is discussed in Sect. 2.2, and the latter in Sect. 2.5. 


Th s large system of equations is solved by a series of 
Newton-Raphson iterations until convergence has been ob- 
tained. This requires the knowledge of the partial derivatives 
of p', r, L' and R! with respect to p.T.L and R. Some of 
these derivatives may be evaluated analytically, but for others 
(such is those involving the equation of state or opacity tables) 
numerical derivatives have to be employed. Each iteration re- 
quires the solution of a linear system of equations, through the 
inversion of a block diagonal matrix. It is important in such a 
relaxation scheme to have initial conditions sufficiently close 
to the true solution; two Runge-Kutta integrations, one outward 
from he centre and one inward from the surface, are used to 
provide such an initial trial solution. 


2.3. Boundary' conditions 

Boundary conditions are imposed at the centre and at the surface. 
At the centre, 
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2.2. Spatial integration 

The equations describing the equilibrium structure of a solar 
model are fourth order, with two boundary conditions each at the 
centre and surface. These were originally solved by the shooting 
technique, whereby the equations were integrated inwards from 
the surface and outwards from the centre, the matching being 
achieved bv varying the four free parameters (central pressure 
and temperature, surface radius and luminosity). This technique 
has largely been superceded by the Henyey relaxation method. 

An iV-point grid m, (i = 1 N with m ; v = M&) is 

imposed, and the differential equations are replaced by a set 
of finite-difference equations. Originally these finite-difference 
equations were second order, but MoSEC uses a fourth-order 
scheme described by Cash & Moore (1980). A - 1 further grid 
points m 3 ,, . . . m v _ 1/2 are interleaved midway between pairs 
of the original X grid points. Equations are then constructed tor 
each of the dependent variables, pressure ip), temperature (2 ). 
luminosity (L) and radius (R)'. for example. 


Pi-t-l - Pi = g ( m J+l “ + + Pi~l) ’ 1 ' ' 

where ' denotes the derivative with respect to the independent 
variable, m. p t + 1/2 * s cH ven 


Similar equations are used for T, L and R. The derivatives 
on the right-hand sides of Eqs (1 ) and (2) are known functions 
of p. T. L and R through the equations of stellar structure. 


wheie £ denotes the energy-generation rate per unit mass, and 
£ eqb the equilibrium energy-generation rate at the centre. At the 

surface, 

Tjv =T(rb.T eff ), (5) 

p iV :=p(Tb,X eff ), (6) 

whe-e T eff is the effective temperature, given by T e 4 ff = 
L v irtoR% (cr is the Stefan-Boltzmann constant), and Tb is 
the optical depth at which interior solution is matched to the 
atmosphere. The latter should be made sufficiently large that 
the diffusion approximation (upon which the radiative-transfer 
equation is based) is valid (Morel et al. 1994). Assuming the 
HSUA X(r, T e ff) law, as fitted by Ando & Osaki (1975). 


r(i.Teff) = -Jett l T + L01T - 0.3 exp( -2.547) 

-0.291 exp(-30r)i , < 7) 

7(-b.r e ff) is obtained directly, and p(iy,.T e {t) is obtained by 
inn grating for the structure of the atmosphere using a Runge- 
Ku ta scheme with adaptive step size. 


2.4. Truncation error 

If h is the step size in mass, then the local truncation error in 
Eq. (1) is 0{h 5 ). This gives a global truncation error which 
scales like h i , so the finite-difference scheme is fourth order. 
This is demonstrated in Fig. 1, which plots the logarithm of the 
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log N 


Fig. 1. Logarithm of the relative errors in radius (solid line) and lu- 
minosity (dashed line) plotted against the logarithm of the number of 
radial grid points, N. 


relative errors in radius and luminosity against log N. The slight 
departure of the errors in luminosity from following a straight 
line is probably due to errors arising from the interpolation of 
opacity and equation of state tables. A value of N = 1000 is 
sufficient to constrain the relative truncation errors in radius and 
luminosity to 10 -6 . 


2.5. Nuclear burning 

The treatment of the time evolution of chemical composition 
is crucial to the success of a stellar evolution code. The princi- 
pal energy-generating cycle of the proton-proton chain (which 
dominates in the sun) consists of two halves, the conversion 
of hydrogen into 3 He via deuterium, and the fusion of 3 He to 
form 4 He. In the centre of the sun, the second of these reactions 
proceeds very much faster than the first, so that equilibrium of 
3 He is quickly established. Away from the centre, equilibrium 
takes much longer to be reached, so it is essential to follow the 
evolution of the 3 He abundance correctly in order to calculate 
the total luminosity accurately. Since the time scales of the two 
reactions of the principal branch of the proton-proton chain are 
so different, the equations constitute a stiff system, and should 
be integrated by an appropriate implicit technique. 

Christensen-Dalsgaard (1991a) used a backward Euler 
method for evolving the 3 He abundance, and a second-order 
scheme for the evolution of the hydrogen abundance. Sub- 
sequent codes have employed more sophisticated integration 
schemes; for instance, the CESAM code developed by Morel 
et al. (1990) uses an implicit, second-order, time-integration 
scheme. MoSEC uses a similar second-order implicit method 
for determining the time evolution of the abundances of H. 3 He, 
4 He, 12 C, 13 C, 14 N, and 16 0. To simplify the reaction network. 
'Be is assumed to be always in equilibrium, while only the 
dominant branches are considered in the CNO cycle. 

The reaction rate r p = r p (m l . t ). corresponding to process 
p (as defined in Table 1) at the point m = m ? , is then given by 


Table 1. Simplified nuclear reaction network. 


p Reaction 


1 1 H(p. /3 + i/) 2 H(p, 7 ) 3 He 

2 3 He( 3 He,2p) 4 He 

3 3 He(Q,7) 7 Be(,3-,i/) 7 Li(p,Q) 4 He 

4 12 C(p, 7 ) 13 N(,^ + ) 13 C 

5 13 C(p,7) l4 N 

6 14 N(p. 7) 15 0(, i'.3 + ) 15 N(p. q) 12 C 

7 16 Q(p,7) 17 F(. i/i 3 + ) 17 Q(p.a) 14 X 


one of the following equations (in units of s *g *): 


T\ = C\X\X\ , 

(8) 

r 2 = c 2 X 3X3 , 

( 9 ) 

r 3 = C3X3X4 , 

(10) 

r 4 = c 4 X\X\2 , 

(ID 

r 5 = C$XiXi 3 , 

(12) 

^6 = ^6^1X14 , 

( 13 ) 

r 7 = cjX\X\q , 

( 14 ) 


where X n = is the fractional abundance by mass of 

the element with atomic number n 7 and c p = c p {rrii . t) is the 
reaction-rate coefficient for process p. The latter are functions of 
density and temperature alone. The evolution of the elemental 
abundances is then calculated according to the following set of 
equations: 

= ai (2r 2 - 3r i - r 3 - r 4 - r 5 - 2r 6 - 2 r 7 ) , (15) 

= a 3 (ri - 2r 2 - r 3 ) , (16) 

= a 4 (r 2 4 - r 3 + r 6 4- r 7 ) , (17) 

= 0-12 (re ~ r 4 ) , (18) 

= at 3 (r 4 -r 5 ), (19) 

= au ( r 5 - r 6 4 - r 7 ) , ( 20 ) 

= ai6(-r 7 ) • (21) 

w here a n is the atomic mass of the element w ith atomic number 
n. 

2.6. Diffusion and gravitational settling 

The relative abundances of chemical elements in the solar in- 
terior are modified not only by nuclear burning, but also by 
the motion of atoms of different species with respect to one 
another. Gravitational and radiative forces acting on individual 
atoms drive such motions, while interactions between atoms 
tend to redistribute momentum in a random way and thereby 


dXi 
d t 
dX 3 
d t 

dX 4 

d t 
dX 12 

dt 

dX l3 

dt 

dXu 

dt 

dAig 

dt 
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counteract these forces. The details of element segregation in 
the sun are determined by the exact nature of the competition 
between these processes in a multi-component plasma. 

Since the abundances of the heavy elements are small, it is a 
reasonable approximation to treat them as if they were trace ele- 
ments in a H-He background (Proffitt 1994). The seven species 
discussed in the context of nuclear reactions (H, 3 He, 4 He, 12 C, 
13 C, 14 N, and 16 0) are treated individually in the settling and 
diffusion calculation, and are assumed to be fully ionized, while 
all the other heavy elements are grouped together and assumed 
to diffuse like fully ionized iron. Since the diffusion of heavy 
elements has a relatively small effect on the properties of the 
computed model, this represents a good approximation (Bahcall 
& Pinsonneault 1995). 

The simplified transport equations of Michaud & Proffitt 
(1993) are used to describe the gravitational settling and diffu- 
sion of helium and the heavy elements. The diffusion velocities 
given by these equations are accurate to within about 10% when 
compared to the more rigorous derivation of Burgers (1969). 
The equations are solved using an explicit, second-order, time 
integration scheme, with time steps determined by the overall 
truncation error requirement. 


2. 7. Time evolution 


An Q-point grid in time, tj , is constructed such that t\ ~ 0 and 
tQ — where is the age of the sun. In order to evolve the 
solar model from time tj to time tj+\, assuming the structure 
Pij = ptj = p(m,.tj), T itj = T{m ly tj) to be 

known at tj, the following procedure is adopted: 


1. The coefficients c p;l j = c p (m t , tj) are evaluated using the 
known values of and T t j. The coefficients c p -ij+ \ are 
initially assumed to be equal to c p;lJ . 

2. Eqs. (15)-(2l) are integrated assuming the coefficients 
c p {m l .t) take the values 


Cp{m t .t) 


t-tj 
tj + 1 — 


'Cp;i,j+1 + 


Ej- n 


(22) 


j + 1 


3. The structure is evaluated at timefj+i, using the new hydro- 
gen abundance so obtained. This yields a new set of values 

for Cp t ,j + i . 

4. Steps 2 and 3 are repeated until satisfactory convergence has 
been obtained. 



log Q 


Fig. 2. Logarithm of the relative errors in radius (solid line) and lu- 
minosity (dashed line) plotted against the logarithm of the number of 
temporal grid points, Q . 


procedure used to update the chemical composition between 
tj and tj+ 1 is second-order, the truncation error in the final 
model at time £ 0 goes like Q~ 2 . In order to exploit this known 
dependence, a series of increasing trial values of Q is used, 
and the final result is obtained by Richardson extrapolation. In 
order to demonstrate this, Fig. 2 shows how the truncation error 
in radius and luminosity varies with Q\ the logarithm of the 
relative error in each of these quantities is plotted against log Q 
for Q = 10, 20, 40, 80. While the relative error is still larger 
than 1CT 5 for Q — 80, Richardson extrapolation enables the 
radius and luminosity to be estimated to an accuracy of better 
than 1 ) -6 using just Q = 10, 20, 40. 

2.8 . Calibration of mixing length and initial hydrogen 
abundance 

The final values at time £ 0 of the luminosity and radius depend 
on the initial hydrogen abundance and the mixing-length pa- 
rametf r, a. These are modified by a series of Newton-Raphson 
iterations until the luminosity and radius agree with L 0 and 
R(? (the present-day luminosity and radius) to within one part 
in 10 6 Sufficiently many grid points are required in the tempo- 
ral and spatial integrations that their relative truncation errors 
are less than 10 -6 . 


The grid points tj are not distributed uniformly in time: in 
the early stages of evolution they are much more closely spaced 
than later on. The distribution is given by the functional form 

tj _ tanh(4.r - 3) - tanh(-3) 
t? tanh(l) - tanh(-3) 

where x — (j - 1)/(Q - 1). The exact function chosen to 
define the grid in time does not affect the rate of convergence 
of the solution as Q is increased, but does affect the actual 
truncation error for a particular value of Q . The form (23) is 
chosen as it gives a particularly low truncation error. Since the 


3. Ph} sical assumptions 

3.1. Si'lar age and luminosity 

Estimates of the age of the oldest meteorites set a low er bound 
on the age of the solar system as a whole, and on the age of 
the su i in particular. Without access to all the meteoritic data, 
Guent ter (1989) recommended a value of 4.49 Gyr. More recent 
studies suggest the value should be nearer 4.6 Gyr (e.g. Bahcall 
et al. '995), and it is this value which is adopted here. At any 
rate, tl e correct value for the solar age remains one of the more 
uncert tin ingredients of solar modelling. 
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The luminosity is also a somewhat uncertain quantity due to 
the difficulty of calibrating satellite radiometers, and also due 
to its inherent long-term variability. The value adopted here is 
the same as that used by Bahcall et al. (1995), namely 3.844 x 
10 33 ergs -1 . 

5.2. Element abundances 

The relative abundances of elements in the solar interior are 
determined by a combination of meteorite analysis and photo- 
spheric line strength measurement. The former gives the initial, 
homogeneous, composition of the primeval solar nebula, while 
the latter reflects the present-day composition of the outer lay- 
ers. Apart from the depleted elements Li and Be, the agreement 
between the two has improved greatly during recent years as 
improvements have been made in the atomic data which affect 
photospheric abundance measurements. This enables us to have 
greater confidence in our knowledge both of the relative abun- 
dances of the heavy elements, and also of Z/X, the ratio of 
heavy-element abundance to hydrogen abundance. Current best 
estimates are given for these quantities by Grevesse & Noels 
(1993), and are adopted in the present work; the recommended 
value for Z/X is 0.0245. 

3.3. Energy generation 

Thermonuclear energy generation is by the three branches of 
the proton-proton chain and by the CNO cycle. The energy pro- 
duction per unit mass is computed using a subroutine written by 
Bahcall (Bahcall & Pinsonneault 1992a, 1992b), with cross sec- 
tions taken from Bahcall & Pinsonneault (1992a), and energy 
releases for each reaction taken from Bahcall & Ulrich (1988). 

3.4 . Equation of state 

There are two broad approaches to the problem of finding the 
thermodynamic properties of a partially-ionized plasma. The 
chemical-picture approach, of which the Saha equation is an 
example, is based on the principle of free-energy minimiza- 
tion. The partition function, Z, of the plasma is assumed to 
be factorizable into a part corresponding to the internal exci- 
tation states of individual particles, and a part corresponding 
to their translational states. The free energy, — ArTTnZ, where 
k is Boltzmann’s constant, is then minimized with respect to 
variations in the occupation numbers which satisfy appropriate 
stoichiometric constraints. The Saha equation is the simplest 
realization of this procedure, but suffers from the disadvantage 
that it predicts unphysical recombination of ions and electrons 
in the solar centre. There have been several attempts to alle- 
viate this problem. In the CEFF equation of state, Eggleton et 
al. (1973) introduced an extra term into the free energy which 
forced complete ionization in the solar centre, although there 
was no physical justification for this addition. Mihalas et al. 
(1988) considered an occupation-probability formalism which 
effectively provided a density-dependent cut-off in the internal 
part of the partition function. This again predicted almost com- 


Table 2. Relative heavy-element abundances. 


Element 

Relative 
mass fraction 

Relative 
number fraction 

C 

0.1906614 

0.2471362 

N 

0.0558489 

0.0620778 

O 

0.5429784 

0.5283680 

Ne 

0.2105114 

0.1624178 


plete ionization in the solar centre. The Mihalas, Hummer and 
Dappen (MHD) equation of state was published in the form of 
tables for a single heavy-element abundance of Z = 0.02, and 
for a heavy-element mixture of carbon, nitrogen, oxygen and 
iron. 

In this study, a new MHD-like equation of state, MHD-E, 
is calculated using a code written by the author, with energy- 
level data kindly provided by W. Dappen. Relativistic effects 
are included in the calculation of the electron free energy (they 
were not in either OPAL or MHD), which will be seen to have a 
significant effect on the computed models. The relative heavy- 
element abundances are set equal to those used in the calculation 
of the OPAL tables, while the total heavy-element abundance, 
Z, may be varied, overcoming the limitation of the original 
MHD tables. 

The other approach to the equation of state is known as the 
physical picture. It does away w ith the concept of atoms, consid- 
ering only fundamental particles such as nuclei and electrons. 
Interactions between particles are taken into account using the 
techniques of many-body theory. The only realization of these 
ideas has been carried out by the OPAL group at Livermore, with 
the publication of preliminary tables in 1994 and subsequently 
of tables with a finer mesh (Rogers et al. 1996). These tables 
are computed using the Grevesse (1993) abundances for carbon, 
nitrogen and oxygen, with the abundances of all the other heavy 
elements being added to the abundance of neon (see Table 2); 
they cover Z = 0.00, Z = 0.02 and Z = 0.04. Various compar- 
isons have been made between the OPAL and MHD equations 
of state (Dappen etal. 1990, Dappen 1992). with the conclusion 
that they are remarkably similar over the bulk of the solar inte- 
rior, with OPAL also predicting almost complete ionization at 
the solar centre; a full explanation for this has yet to be found. 

Evolutionary models are computed using both equation of 
state formalisms. The pressure and other thermodynamic quan- 
tities are evaluated as functions of the density, temperature, hy- 
drogen and overall heavy-element mass fractions by means of 
interpolating tables. Some studies have ignored the variation of 
the overall heavy-element mass fraction in the calculation of 
the equation of state, using instead a constant, prescribed value 
(e.g. Morel et al. 1997). We perform a comparison of models 
calculated with and without this assumption in order to test its 
significance. In our best model, we choose only to ignore the 
variation of the individual heavy-element abundances. 
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Table 3. Comparison of the properties of solar models. 


Model 

2Cjnit 

a 

Xc 

Fsurf 

Zsnrt 

T c 

(10 6 K) 

Pc 

(gem -3 ) 

$3701 

(SNU) 

$-lGa 

(SNU) 

xcz 

opal 

0.7192 

1.9351 

0.3618 

0.2634 

0.0176 

15.468 

150.79 

6.36 

122.4 

0.726 

ydiff 

0.7173 

2.0983 

0.3463 

0.2365 

0.0183 

15.588 

154.82 

7.22 

127.0 

0.711 

zdiff 

0.7115 

2.0792 

0.3395 

0.2422 

0.0181 

15.671 

154.83 

7.86 

130.4 

0.712 

zdiffop 

0.7116 

2.0754 

0.3396 

0.2421 

0.0181 

15.670 

154.78 

7.84 

130.4 

0.712 

mhd 

0.7092 

2.0383 

0.3380 

0.2439 

0.0181 

15.677 

154.96 

7.93 

130.8 

0.712 

mix 

0.7127 

2.0707 

0.3409 

0.2438 

0.0181 

15.656 

154.71 

7.73 

129.8 

0.712 

zdiff re 

0.7104 

2.0683 

0.3389 

0.2430 

0.0181 

15.672 

154.90 

7.87 

130.5 

0.712 

Proffitt (1994) a 
Bahcall & 

0.6984 

1.711 

0.3288 

0.7290 

0.0196 

15.81 

155.9 

9.02 

136.9 

0.712 

Pinssoneault (1995) b 

0.7025 

2.09 

0.3333 

0.7351 

0.0180 

15.84 

156.2 

9.3 

137 

0.712 

Richard et al. (1996) c 
Christensen-Dalsgaard 

0.7012 

1.768 

0.3328 

0.7226 

0.0190 

15.67 

154.53 

8.49 

132.8 

0.716 

et al. (1996) d 

0.7091 

1.9905 

0.3353 

0.7373 

0.0181 

15.668 

154.24 

8.2 

132 

0.712 

Morel et al. (1997) e 

0.7064 

1.92 

0.3412 

0.737 

0.0180 

15.65 

151.2 

8.27 

130 

0.711 


a CEFF equation of state, age of 4.6 Gyr, Z/X = 0.0269. 
b CEFF equation of state, age of 4.57 Gyr. 
c MHD equation of state, includes rotation induced mixing. 
d OPAL equation of state, age of 4.6 Gyr. 
e OPAL equation of state, age of 4.55 Gyr. 


3.5 . Opacity 

As was the case for the equation of state, the opacity is obtained 
from interpolating tables. The tables are constructed using the 
OPAL opacities (Rogers & Iglesias 1992), calculated with the 
Grevesse (1993) heavy -element mixture, except at low temper- 
atures (;$ 10 4 K), where the Kurucz (1991) low-temperature 
opacities are substituted. The interpolation is carried out using a 
package written by G. Houdek (Houdek & Rogl 1996), which al- 
lows a choice between a minimum-norm and a birational-splines 
algorithm. There is also a freedom of choice in the number of 
fitting points used in the X and Z interpolations (2,3 or 4). In 
the first instance the minimum-norm algorithm is used with four 
fitting points; the effect of other choices on the calibrated hy- 
drogen abundance and mixing-length parameter is investigated 
in Sect. 4.2. 

The heavy-element abundance used to interpolate the opac- 
ity tables only includes the abundance changes due to element 
segregation, since nuclear burning by the CNO cycle has little 
effect on the opacity (Proffitt 1994). As with the equation of 
state, the variation of the individual heavy-element abundances 
is ignored; as suggested by Morel et al. (1997), this may well 
have a significant effect on the computed models relative to the 
high precision of p mode frequency data. 

4. Results 

4. 1. Models constructed 

The MoSEC code is used to construct a series of calibrated, 
evolutionary solar models. For all the models, the matching op- 


tical cepth, Tfe, is chosen to be 2, while a global truncation error 
of bet .er than 10 -6 is sought. From the considerations of Sect. 
2, 1000 radial grid points is sufficient to meet this requirement 
in the space dimension, while Richardson extrapolation with 
Q = 10, 20, 40 provides sufficient accuracy in time. The mod- 
els constructed are as follows: 

- oj al : solar model constructed using the OPAL equation 
of state with no settling or diffusion of helium or heavy 
el ;ments. 

- yc iff : solar model including the settling and diffusion of 
helium, and using the OPAL equation of state. 

- zciff : solar model including the settling and diffusion of 
helium and heavy elements, and using the OPAL equation 
of state. 

- zdiffop : solar model including the settling and diffusion ot 
hi lium and heavy elements, using the OPAL equation of 
st ite w ithout taking into account the variation of the heavy- 
el iment abundance. 

- mad : solar model including the settling and diffusion ot 
h( lium and heavy elements, using the MHD equation ot 
st ite. 

- mix : solar model including the settling and diffusion of he- 
lii im and heavy elements, using the OPAL equation of state, 
at d including a simple model of turbulent mixing below the 
bi .se of the convection zone. 

- zc iffre : solar model including the settling and diffusion 
ol helium and heavy elements, using the OPAL equation of 
st ite, and including the relativistic correction to the electron 
pi essure. 
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Table 4. The effect on the calibrated hydrogen abundance and mixing- 
length parameter of the choice of opacity interpolation method and 
number of interpolation points. 


Interpolation 

method 

Number of 
points 

Ainit 

Q 

minimum norm 

2 

0.718695 

1.94696 


3 

0.719119 

1.94486 


4 

0.719232 

1.94403 

birational splines 

4 

0.719177 

1.94298 


In all cases, the heavy elements are initially chosen to be 
in the Grevesse & Noels (1993) proportions and iteration is 
performed to ensure that the final surface heavy-element mass 
fraction satisfies ZjX — 0.0245 (see Sect. 3.2). It would be pos- 
sible in the case of the models with heavy-element settling and 
diffusion to iterate the individual heavy-element abundances to 
ensure that the final surface abundances are in the Grevesse & 
Noels (1993) proportions, a procedure carried out by Richard et 
al. (1996). However, bearing in mind that neither the equation 
of state nor opacity calculation takes into account the variation 
of individual heavy-element abundances, we do not believe that 
this procedure would have a significant effect, and therefore do 
not include it at present. 

Various properties of the computed models are listed in Ta- 
ble 3 and compared with the properties of a selection of other 
recent solar model calculations which include settling and dif- 
fusion of helium and heavy elements. X c , T c and p c denote the 
central hydrogen mass fraction, central temperature and central 
density respectively. A'i n i t denotes the initial (uniform) hydro- 
gen mass fraction, y sur f denotes the final surface helium mass 
fraction, Z sur f denotes the final surface heavy-element mass 
fraction, while xcz denotes the fractional radius of the base of 
the convection zone. 4 >j7 C i and 4>riG a represent the predicted 
neutrino capture rates for 37 C1 and 71 Ga detectors, given in 
solar neutrino units (1 SNU = 10“ 36 capture atom' 1 s“ l ). 

4.2. Opacity ■ and equation of state interpolation 

A series of models is constructed using the two different opacity 
interpolation methods (minimum norm, and birational splines) 
and various different numbers of fitting points (2,3 and 4); all 
other aspects of the physics (OPAL equation of state, no set- 
tling or diffusion) are kept constant. The calibrated hydrogen 
abundances and mixing-length parameters for these models are 
listed in Table 4. The third model listed corresponds to the “opal" 
model in table 3. 

As can be seen from this table, there is a reasonably good 
agreement between models constructed using the two differ- 
ent interpolation algorithms with four interpolation points. The 
agreement is better for the initial hydrogen abundance than for 
the mixing-length parameter, which is as expected given that 
the latter is more sensitive to the opacity. The three models 


Table 5. The effect on the calibrated hydrogen abundance and mixing- 
length parameter of changing the fineness of the equation of state grid. 


EOS grid 

Ajnit 

a 

coarse 

0.716636 

1.890395 

fine 

0.716632 

1.889529 


constructed with the birational-spline algorithm show, as ex- 
pected, that as the number of interpolation points is reduced, 
the calibrated hydrogen abundance and mixing-length parame- 
ter exhibit increasingly large errors. 

The effect of changing the fineness of the grid used in the 
equation of state tables is show n in Table 5. Both the models 
listed in this table are constructed without helium settling or 
diffusion, using the MHD equation of state; the first is con- 
structed using equation of state tables on the same grid as the 
OPAL equation of state, while the second is a similar model con- 
structed using MHD equation of state tables with a grid twice 
as fine. The differences between the two models in this case are 
smaller (much smaller in the case of the hydrogen abundance) 
than the differences described above between models calculated 
using two different opacity interpolation schemes. The errors in- 
troduced by interpolation (especially in the calibrated hydrogen 
abundance) are therefore dominated by the those arising from 
the opacity interpolation. The original choice of grid for the 
MHD equation of state tables (corresponding to that used in the 
OPAL tables) is therefore sufficiently fine. 

4.5. Microscopic diffusion 

It is well known that models without settling and diffusion of 
helium do not very well reproduce the seismically determined 
surface helium abundance and convection-zone depth of the 
sun (e.g. Christensen-Dalsgaard et al. 1993, Proffitt 1994), hav- 
ing surface helium abundances too high, and convection-zone 
depths too low. In this study, including settling and diffusion of 
helium decreases Y su r f from 0.2634, in the case of the “opaP 
model, to 0.2365, in the case of the “ydiff" model, and decreases 
xcz from 0.726 to 0.71 1. These are similar findings to those of 
Christensen-Dalsgaard et al. (1993). Comparison may be made 
writh the seismically determined values of ~ 0.25 (Diippen et al. 
1988, Vorontsov et al. 1991, Basu & Antia 1995, Dziembowski 
et al. 1994) for >' sur f, and 0.713±0.003 (Christensen-Dalsgaard 
et al. 1991b) for xcz- 

In addition to improving the agreement in convection-zone 
depth and surface helium abundance, helium settling and diffu- 
sion dramatically improve the agreement in sound-speed with 
the sun. This is shown in Fig. 3, where the dashed line represents 
the relative difference in squared sound speed. Sc 2 /c 2 . between 
the “opal" model and the sun, and the dotted line represents 
the same quantity for the “ydiff" model. This improvement is 
principally due to the increased convection-zone depth, but also 
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Fig. 3. The difference in squared sound speed between the models zdiff, 
ydiff and opal and the sun, plotted as a function of fractional solar 
radius, constructed from an inversion of 2 months of MDI p mode 
frequency data with model S of Christensen-Dals gaard et al. (1996) as 
reference. The thin solid lines show the standard errors in the sound- 
speed inversion. 

partly due to the reduced mean molecular weight just below the 
base of the convection zone. There are still significant features 
just below the base of the convection zone, where the sound 
speed is too low in the model, and near the centre, where the 
sound speed is too high in the model. 

Heavy-element settling and diffusion have been included in 
many of the recent solar model calculations, including all the 
models listed in the lower half of Table 3. These processes are 
included in model “zdiff’, increasing V; urf by 0.0057 and re- 
ducing xqz by 0.001, in very close agreement with the results 
of Proffitt (1994). The corresponding form of <5c 2 /c 2 , showm by 
the thick solid line in Fig. 3, is very similar to that found for 
the “ydiff’ model, with a maximum difference of about 0.05% 
in sound speed between the two models at around r/i£ 0 — 0.5 
(which is, nevertheless, significant compared to the errors in 
the inversion, shown by the thin solid lines). Consequently, un- 
like helium, heavy-element settling and diffusion do not sig- 
nificantly improve the agreement in sound speed between solar 
models and the sun. 

The flux of neutrinos arising from the decay of 8 B is strongly 
affected by the inclusion in the evolution calculation of the dif- 
fusion and settling of helium and heavy -elements. Gravitational 
settling leads to an increase in mean molecular weight at the cen- 
tre. and a consequent increase in the central temperature. This 
leads to an increase in the flux of S B neutrinos (since the latter 
goes roughly like I 16 at the solar centre). Bahcall & Pinson- 
neault (1992b) found that helium settling alone brought about 
an increase of about 12% in the flux measured by the chlorine 
neutrino capture experiments: Proffitt (1994) found an increase 
in this flux of about 13.6%, the discrepancy being due to the 
different formulations of the diffusion equations for helium. In 
this study, we find that, the 37 C1 flux goes from 6.12 SNU to 
6.95 SNU, an increase of about 14%: this is very close to the 
value obtained by Proffitt (1994). 


Fig. 4. The relative difference in squared sound speed between the 
models zdiff, zdiffop, mix, mhd and zdiffre and the sun, constructed in 
the same way as in Fig. 3. The thin solid lines represent the standard 
errors in the sound-speed inversion. 

We find that helium and heavy-element settling combined 
produce an increase in the 3/ Cl flux of about 24%. There is 
considt rable variation in the corresponding increase found by 
other authors, varying from 15% (Richard et al. 1996) to 32% 
(Bahcall & Pinsonneault 1995); Turck-Chieze & Brun (1997) 
find an increase of 26%, closest to our result. The reasons for 
these wide discrepancies are not yet clear. In terms of the actual 
37 C1 flux prediction, all the models listed in Table 3 predict 
a somewhat higher 37 C1 flux than the model “zdiff’, although 
the calculations of Turck-Chieze & Brun (1997) predict a lower 
37 C1 flux of 7.2 SNU. 

4.4. The equation of state 

The model “zdiffop” is constructed only taking into account 
the variation of the overall heavy-element abundance, Z, in the 
interpolation of the opacity tables. This is a similar assumption 
to that made by Morel et al. (1997). The relative difference in 
squarei l sound speed between this model and the sun is shown by 
the dot ed line in Fig. 4, and compared with the same quantity for 
four ot uer models, “zdiff’ (solid line), “mhd” (dot-dashed line), 
“mix” i dashed line) and “zdiffre” (dot-dot-dot-dashed line). The 
relative differences in squared sound speed between the models 
“zdiffop”, “mix”, “mhd” and “zdiffre” and the model “zdiff* 
are plotted in Fig. 5 using the same line styles. As can be seen, 
the difference in Sc 2 /c 2 between “zdiffop” and “zdiff’ is almost 
within the standard errors in the inversion (shown by the thin 
solid li les), and the neglect of the variation of the overall heavy- 
elemert abundance in the calculation of the equation of state is 
therefc re valid. It therefore seems that, contrary to the claims of 
Morel st al. (1997), equation of state data taking into account 
detaile I changes in the heavy-element mixture are unnecessary. 

Conparison is made between the OPAL and the MHD-E 
equations of state using the “mhd” model, which differs from 
the “zdiff’ model only in using the MHD-E instead of the OPAL 
equation of state. Early comparisons of the OPAL and MHD 
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Fig. 5. The relative difference in squared sound speed between the 
models zdiffop, mix, mhd and zdiffre and the model zdiff. The thin 
solid lines represent the standard errors in the sound-speed inversion. 

equations of state (Dappenetal. 1990, Dappen 1992) found very 
close agreement over most of the solar interior. More recently, 
Guenther et al. (1996) computed models using both equation of 
state formalisms, finding very close agreement both in global 
properties and in local physical quantities, while Rogers et al. 
(1996) found slight differences in the location of the hydro- 
gen and helium ionization zones leading to a so-called ‘"pres- 
sure spike” in the MHD equation of state at temperatures of 
~ 50 000 K. 

In the present calculation, the two models “mhd” and “zd- 
iff” do not agree as closely as those of Guenther et al. (1996), 
with larger differences in initial hydrogen abundance, mixing- 
length parameter, and also in the sound speed. However, better 
agreement is found when the “mhd” model is compared to the 
“zdiffre” model, which uses the OPAL equation of state but in- 
cludes a correction for relativistic effects. This shows that the 
inclusion of the relativistic correction in the new MHD-E equa- 
tion of state accounts for much of the difference found between 
the “mhd” and “zdiff” models. The relativistic correction has 
the effect of increasing the electron pressure at the centre, thus 
reducing the hydrogen abundance there in the calibrated model. 
This has a significant effect on the sound speed throughout the 
radiative interior, reducing it by a maximum of about 0.1 % near 
the base of the convection zone. 

4.5. Sub convection-zone mixing 

A simple model of the sub convection-zone mixing associated 
with the solar tachocline (Spiegel & Zahn 1992, Elliott 1997) 
is included in the model “mix”. Instead of assuming complete 
mixing to the base of the convection zone, as in the standard solar 
model prescription, the region of mixing is assumed to extend a 
distance of O.O 2 i ?0 below the base of the convection zone: this 
value is somew hat low'er than the thickness of the tachocline as 
inferred from rotation-rate inversions (Kosovichev 1996). The 
effect of the extra mixing is to redistribute gravitationally set- 
tled helium back into the convection zone, thereby reducing the 


mean molecular weight just below the base of the convection 
zone. This, in turn, increases the sound speed in that region as 
is shown by the dashed line in Fig. 5. Deeper down, the recal- 
ibration of the model reverses the situation, causing the sound 
speed to be lower in “mix” than in “zdiff’. The overall effect of 
the extra mixing is to improve slightly the agreement in sound 
speed with the sun: the sound speed increase in the “mix” model 
relative to the “zdiff’ model would need to be displaced down- 
ward to account fully for the hump in 5c 2 /c 2 just below the base 
of the convection zone. 


5. Conclusion 

The new solar evolution code described in this paper has been 
shown to produce results broadly consistent with previous cal- 
culations, as well as having well-controlled truncation errors 
both in space and time. The effects of using interpolating tables 
for the opacity and equation of state have been quantified, with 
the conclusion that at present the interpolation of the OPAL 
opacity tables sets the limit on errors in the calibrated hydrogen 
abundance and mixing-length parameters. 

The effect of including helium diffusion and settling in the 
evolution calculation is very similar to that found by previous 
authors. There is a significant increase in the predicted neu- 
trino flux for the chlorine experiments (14%), and a somewhat 
smaller increase (3.8%) for the gallium experiments. When 
heavy-element diffusion is also included, the corresponding in- 
creases relative to the non-diffusive models are 23% and 6 . 5 %, 
both similar to the values found by Turck-Chieze & Brun ( 1 997). 
The actual value of the predicted count rate for the chlorine ex- 
periment, 7.86 SNU, is somewhat smaller than that reported by 
most other authors, but larger than the 7.2 SNU reported by 
Turck-Chieze & Brun. These differences point to the need for 
further investigation. 

The inclusion of the variation of the overall heavy-element 
mass fraction in the calculation of the equation of state has been 
shown to have an almost negligible effect on the sound speed 
on the computed model, rendering the effect of the variation of 
individual heavy-element abundances in the equation of state 
calculation undetectable at the current precision of structure 
inversions. This is, however, not the case in the calculation of the 
opacity, and opacity data taking into account detailed changes 
in composition are necessary. 

A comparison is provided here between models constructed 
using the OPAL and MHD-E equations of state. The agreement 
is not as good as that reported in previous comparisons between 
OPAL and MHD, principally due to inclusion of the relativis- 
tic correction to the electron pressure in the new'Iv-calculated 
MHD-E equation of state. The significant effect that this correc- 
tion has on the sound speed of calibrated models should perhaps 
prompt efforts to include it consistently in future equation of 
state calculations. 

A simple model of the sub convection-zone mixing associ- 
ated with the solar tachocline has been shown to improve the 
agreement in sound speed with the sun just below the base of the 
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convection zone. Further study is required to understand fully 
the sound speed enhancement found in this region of the sun. 

In spite of improvements in solar evolution calculations, 
it is as well to point out that there are still considerable un- 
knowns in many of the physical assumptions (Gough & Toomre 
1991). The high accuracy of the calculation presented here fa- 
cilitates detailed comparison with models by other authors, but 
belies some of the inherent uncertainties. Perhaps most serious 
is the assumption of uniform chemical composition at zero age, 
which is contingent on the sun having passed through a fully- 
convective Hayashi phase of evolution. Calculations by Larson 
(1969) have suggested that this may not have occurred. Another 
significant failing of the standard solar model is its neglect of 
large-scale macroscopic motions (aside from convection). It is 
known that such motions must occur in any rotating star, which 
could have important repercussions for its evolution by redis- 
tributing both material and angular momentum. Mass loss and 
accretion are known to be important processes in stellar evolu- 
tion but are again entirely neglected in the standard solar model 
prescription. The modular nature of MoSEC should permit re- 
finements of the physical assumptions to be easily incorporated, 
allowing the code to keep pace with our improving understand- 
ing of solar evolution. 
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